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Abstract 

We give an overview of recent developments in the problem of reconstructing 
a band-limited signal from non-uniform sampling from a numerical analysis view 
point. It is shown that the appropriate design of the finite-dimensional model plays 
a key role in the numerical solution of the non-uniform sampling problem. In the one 
approach (often proposed in the literature) the finite-dimensional model leads to an 
ill-posed problem even in very simple situations. The other approach that we con- 
sider leads to a well-posed problem that preserves important structural properties 
of the original infinite-dimensional problem and gives rise to efficient numerical al- 
gorithms. Furthermore a fast multilevel algorithm is presented that can reconstruct 
signals of unknown bandwidth from noisy non- uniformly spaced samples. We also 
discuss the design of efficient regularization methods for ill-conditioned reconstruc- 
tion problems. Numerical examples from spectroscopy and exploration geophysics 
demonstrate the performance of the proposed methods. 
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1 Introduction 



The problem of reconstructing a signal / from non-uniformly spaced measurements f(tj) 
arises in areas as diverse as geophysics, medical imaging, communication engineering, 
and astronomy. A successful reconstruction of / from its samples f(tj) requires a priori 
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information about the signal, otherwise the reconstruction problem is ill-posed. This a 
priori information can often be obtained from physical properties of the process generating 
the signal. In many of the aforementioned applications the signal can be assumed to be 
(essentially) band-limited. 

Recall that a signal (function) is band-limited with bandwidth Q if it belongs to the 
space B n , given by 

B n = {fe L 2 (R) : f(u) = for \u\ > ft} , (1) 

where / is the Fourier transform of / defined by 

+00 

/>)= f f{t)e~ M dt. 



For convenience and without loss of generality we restrict our attention to the case ft — |, 
since any other bandwidth can be reduced to this case by a simple dilation. Therefore we 
will henceforth use the symbol B for the space of band-limited signals. 

It is now more than 50 years ago that Shannon published his celebrated sampling 
theorem ||35|| . His theorem implies that any signal / E B can be reconstructed from its 



regularly spaced samples {/(n)} ne z by 

/W = £ /( „)!EI&Z2>. (2) 

In practice however we seldom enjoy the luxury of equally spaced samples. The solu- 
tion of the nonuniform sampling problem poses much more difficulties, the crucial ques- 
tions being: 

• Under which conditions is a signal f E B uniquely defined by its samples {f(tj)}jez? 

• How can / be stably reconstructed from its samples f(tj)7 

These questions have led to a vast literature on nonuniform sampling theory with deep 
mathematical contributions see |TT], |25|, |3|, |6|, [H| to mention only a few. There is also no 



lack of methods claiming to efficiently reconstruct a function from its samples 0, f|l], [l] 



H|, pED|, |2T| [THj . These numerical methods naturally have to operate in a finite-dimensional 



model, whereas theoretical results are usually derived for the infinite-dimensional space 
B. From a numerical point of view the "reconstruction" of a bandlimited signal / from 
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a finite number of samples {/(^)K=i amounts to computing an approximation to / (or 

/) at sufficiently dense (regularly) spaced grid points in an interval (ti,t r ). 

Hence in order to obtain a "complete" solution of the sampling problem following 
questions have to be answered: 

• Does the approximation computed within the finite-dimensional model actually con- 
verge to the original signal /, when the dimension of the model approaches infinity? 

• Does the finite-dimensional model give rise to fast and stable numerical algorithms? 

These are the questions that we have in mind, when presenting an overview on recent 
advances and new results in the nonuniform sampling problem from a numerical analysis 
view point. 

In Section || it is demonstrated that the celebrated frame approach does only lead to 
fast and stable numerical methods when the finite-dimensional model is carefully designed. 
The approach usually proposed in the literature leads to an ill-posed problem even in very 
simple situations. We discuss several methods to stabilize the reconstruction algorithm 
in this case. In Section |3] we derive an alternative finite-dimensional model, based on 
trigonometric polynomials. This approach leads to a well-posed problem that preserves 
important structural properties of the original infinite-dimensional problem and gives rise 
to efficient numerical algorithms. Section ^ describes how this approach can be modified in 
order to reconstruct band-limited signals for the in practice very important case when the 
bandwidth of the signal is not known. Furthermore we present regularization techniques 
for ill-conditioned sampling problems. Finally Section [| contains numerical experiments 
from spectroscopy and geophysics. 

Before we proceed we introduce some notation that will be used throughout the paper. 
If not otherwise mentioned \\h\\ always denotes the L 2 (R)-norm (£ 2 (Z)-norm) of a function 
(vector). For operators (matrices) ||T|| is the standard operator (matrix) norm. The 
condition number of an invertible operator T is defined by k(A) = \\A\\ and the 
spectrum of T is cr(T). I denotes the identity operator. 



1.1 Nonuniform sampling, frames, and numerical algorithms 

The concept of frames is an excellent tool to study nonuniform sampling problems [13[ 



[I], [15], ^3 ] . The frame approach has the advantage that it gives rise to deep theoretical 
results and also to the construction of efficient numerical algorithms - if (and this point 
is often ignored in the literature) the finite-dimensional model is properly designed. 
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Following Duffin and Schaeffer [n|, a family {ft}j^z i n a separable Hilbert space H is 
said to be a frame for H, if there exist constants (the frame bounds) A, B > such that 

^ll/ll 2 <EK«>l 2 ^ll/ll 2 ' vfeH. (3) 

3 

We define the analysis operator T by 

T:feH^Ff = {(f,fj)} j& , (4) 
and the synthesis operator, which is just the adjoint operator of T, by 

T*:ce£ 2 (Z)^T*c = Y,c 3 f J . (5) 

3 

The frame operator S is defined by S = T*T, hence Sf = ^2j(f, ft) ft- S is bounded by 
AI < S < BI and hence invertible on H. 

We will also make use of the operator TT* in form of its Gram matrix representation 
R : £ 2 (Z) -> £ 2 (Z) with entries R u = (ft, /,). On ft(T) = 72.(i2) the matrix J2 is bounded 
by AI < R < BI and invertible. On £ 2 (Z) this inverse extends to the Moore-Penrose 
inverse or pseudo-inverse R + (cf. ||12||). 

Given a frame {/jjjez for ff, any f E H can be expressed as 

where the elements 7^ := S~ x fj form the so-called dual frame and the frame operator 
induced by jj coincides with S 1-1 . Hence if a set {fj}j^z establishes a frame for H, we 
can reconstruct any function / E H from its moments (/, ft). 

One possibility to connect sampling theory to frame theory is by means of the sine- 
function 

, s sin nt 

smc t = — - . 7 
nt 

Its translates give rise to a reproducing kernel for B via 

f(t) = {f,smc(--t)) VtjEB. (8) 
Combining (H) with formulas (0) and @ we obtain following well-known result [13 
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Theorem 1.1 If the set {sinc(- — tj)}j £ z is a frame for B, then the function f G B is 
uniquely defined by the sampling set {f{tj)}j^z- In this case we can recover f from its 
samples by 

fit) = ^ fihhj , where ^ = S^sinc^ - tj) , (9) 

or equivalently by 

f(t) = CjSincit — tj) , where Rc = b , (10) 

with R being the frame Gram matrix with entries Rji = sinc(tj — t{) and b = {bj} = 
{.A/,))- 

The challenge is now to find easy-to-verify conditions for the sampling points tj such 
that {sinc(- — £j)}jez (or equivalently the exponential system {e 2mtjtJ }j^z) is a frame for 
B. This is a well-traversed area (at least for one-dimensional signals), and the reader 
should consult |], |24| for further details and references. If not otherwise mentioned 



from now on we will assume that {sinc(- — £j)}iez is a frame for B. 

Of course, neither of the formulas (|9|) and (0) can be actually implemented on a 
computer, because both involve the solution of an infinite-dimensional operator equation, 
whereas in practice we can only compute a finite-dimensional approximation. Although 
the design of a valid finite-dimensional model poses severe mathematical challenges, this 
step is often neglected in theoretical but also in numerical treatments of the nonuniform 
sampling problem. We will see in the sequel that the way we design our finite-dimensional 
model is crucial for the stability and efficiency of the resulting numerical reconstruction 
algorithms. 

In the next two sections we describe two different approaches for obtaining finite- 
dimensional approximations to the formulas @ and fllPl). The first and more traditional 
approach, discussed in Section |2], applies a finite section method to equation (p^0|). This 
approach leads to an ill-posed problem involving the solution of a large unstructured linear 
system of equations. The second approach, outlined in Section |3|, constructs a finite model 
for the operator equation in @ by means of trigonometric polynomials. This technique 
leads to a well-posed problem that is tied to efficient numerical algorithms. 

2 Truncated frames lead to ill-posed problems 

According to equation ( JTOj ) we can reconstruct / from its sampling values f(tj) via fit) = 
Yljez c j sinc(t — tj), where c = R + b with bj = f(tj),j G Z. In order to compute a finite- 
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dimensional approximation to c = {cj}j E % we use the finite section method JT7]. For 
x E £ 2 (Z) and n E N we define the orthogonal projection P n by 

PnX (. . . , 0, 0, X— n , X_ n _|_i, . . . , X n —i, x n , 0, 0, ... ) (H) 

and identify the image of P n with the space C 2n+1 . Setting R n = P n RP n and = P n b, 
we obtain the n-th approximation to c by solving 

R nC ^) = &(») . (12) 

It is clear that using the truncated frame {sinc(- — t,)}™ = _ n in C[T0|) for an approximate 
reconstruction of / leads to the same system of equations. 

If {sinc(- — tj)}j e z is an exact frame (i.e., a Riesz basis) for B then we have following 
well-known result. 

Lemma 2.1 Let {sinc(- — tj)}j e % be an exact frame for B with frame bounds A,B and 

Rc = b and R n c^ = b^ as defined above. Then i?" 1 converges strongly to and hence 
c O) _> c f or n _>. qo _ 

Since the proof of this result given in M is somewhat lengthy we include a rather short 
proof here. 

Proof: Note that R is invertible on £ 2 (Z) and A < R < B. Let x E C 2n+1 with ||a;|| = 1, 
then (R n x,x) = (P n RP n x,x) = (Rx,x) > A. In the same way we get \\R n \\ < B, hence 
the matrices R n are invertible and uniformly bounded by A < R n < B and 

4 < R' 1 < 4 for all n E N. 
B ~ A 



The Lemma of Kantorovich [32] yields that R n 1 — > R 1 strongly. □ 



If {sinc(- — tj)}j(zz is a non-exact frame for B the situation is more delicate. Let us 
consider following situation. 

Example 1: Let f E B and let the sampling points be given by tj = ^, j E Z, 1 < m E N, 
i.e., the signal is regularly oversampled at m times the Nyquist rate. In this case the 
reconstruction of / is trivial, since the set {sinc(- — tj)}j^z is a tight frame with frame 
bounds A = B = m. Shannon's Sampling Theorem implies that / can be expressed as 
f(t) = X]jez c i smc (^ — A?) wnere c j — an d the numerical approximation is obtained 
by truncating the summation, i.e., 



m 

j=-n 



fn(t)= V l^-sincit-tj, 

' ' rn 
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Using the truncated frame approach one finds that R is a Toeplitz matrix with entries 



sin — (j — I) 

in other words, i? n coincides with the prolate matrix |36]. [39|. The unpleasant numerical 
properties of the prolate matrix are well-documented. In particular we know that the sin- 
gular values X n of R n cluster around and 1 with log n singular values in the transition re- 
gion. Since the singular values of R n decay exponentially to zero the finite-dimensional re- 
construction problem has become severely ill-posed |P2| |, although the infinite-dimensional 



problem is "perfectly posed" since the frame operator satisfies S = ml, where I is the 
identity operator. 

Of course the situation does not improve when we consider non-uniformly spaced 
samples. In this case it follows from standard linear algebra that o~(R) C {0 U [A, B]}, or 
expressed in words, the singular values of R are bounded away from zero. However for 
the truncated matrices R n we have 

a(i? n )C{(0,£?]} 

and the smallest of the singular values of R n will go to zero for n — > oo, see |[23|| . 

Let A = UT,V* be the singular value decomposition of a matrix A with X = diag({Afc}). 
Then the Moore-Penrose inverse of A is A + = VH + U* , where (e.g., see fll8| ) 

E- = diag({AJ}), K=\T k if t| W °' (13) 

I (J otherwise. 

For R n = U n Tj n V n this means that the singular values close to zero will give rise to 
extremely large coefficients in R+. In fact > oo for n — > oo and consequently 

does not converge to c. 

Practically \\Rn\\ is always bounded due to finite precision arithmetics, but it is clear 
that it will lead to meaningless results for large n. If the sampling values are perturbed due 
to round-off error or data error, then those error components which correspond to small 
singular values are amplified by the (then large) factors l/A^. Although for a given 
R n these amplifications are theoretically bounded, they may be practically unacceptable 
large. 

Such phenomena are well-known in regularization theory fl2| . A standard technique 



to compute a stable solution for an ill-conditioned system is to use a truncated singular 
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value decomposition (TSVD) JT2]. This means in our case we compute a regularized 
pseudo-inverse R^' T = K l E+' r £7* where 

E- = diag(K}), 4=ll /Xk [ \t"~ T ' < 14 > 

otherwise. 



In |23| it is shown that for each n we can choose an appropriate truncation level r such 
that the regularized inverses R^ ,T converge strongly to R + for n — > oo and consequently 
lim ||/ — /' n )|| = 0, where 

n^oo 

n 

/(»>(*) = ^ cf r) sinc(t - t,) 

with 

The optimal truncation level r depends on the dimension n, the sampling geometry, and 
the noise level. Thus it is not known a priori and has in principle to be determined for 
each n independently. 

Since r is of vital importance for the quality of the reconstruction, but no theoretical 
explanations for the choice of r are given in the sampling literature, we briefly discuss 
this issue. For this purpose we need some results from regularization theory. 

2.1 Estimation of regularization parameter 

Let Ax = y 5 be given where A is ill-conditioned or singular and y s is a perturbed right- 
hand side with \\y — y s \\ < 5\\y\\. Since in our sampling problem the matrix under consid- 
eration is symmetric, we assume for convenience that A is symmetric. From a numerical 
point of view ill-conditioned systems behave like singular systems and additional infor- 
mation is needed to obtain a satisfactory solution to Ax = y. This information is usually 
stated in terms of "smoothness" of the solution x. A standard approach to qualitatively 
describe smoothness of x is to require that x can be represented in the form x = Sz with 
some vector z of reasonable norm, and a "smoothing" matrix S, cf. [0, ^9[. Often it is 
useful to construct S directly from A by setting 

S = A p , peN . (15) 
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Usually, p is assumed to be fixed, typically at p = 1 or p = 2. 

We compute a regularized solution to Ax = y 5 via a truncated SVD and want to 
determine the optimal regularization parameter (i.e., truncation level) r. 

Under the assumption that 

x = Sz, \\Ax -y 5 \\ < A\\z\\ (16) 



it follows from Theorem 4.1 in [29] that the optimal regularization parameter r for the 
TSVD is 

*=r^y\ (it) 



where 7i = 72 = 1 (see Section 6 in p9[ ). 

However z and A are in general not known. Using ||Ar— < and ||y|| = \ Ax\ = 

\\ASz\\ = \\A p+1 z\\ we obtain \\y\\ < ||A|| P+1 ||2;||. Furthermore, setting 5\\y\\ = A\\z\ 
implies 



A<5\\A\\ P+1 
Hence combining (|1~?D and (p~8|) we get 



i <»*Hwf*)* (19) 



P J \p 

Applying these results to solving R n c^ = via TSVD as described in the previous 
section, we get 

1 1 1 

S \ P+ 1 /5\ p+1 (8 \ p+1 



f<\M[~) <\\R\\{-) =B\^) , (20) 

where B is the upper frame bound. Fortunately estimates for the upper frame bound are 
much easier to obtain than estimates for the lower frame bound. 

Thus using the standard setting p=lorp = 2a good choice for the regularization 
parameter r is 

tC[B{S/2)^ 3 ,B(S)^ 2 ]. (21) 
Extensive numerical simulations confirm this choice, see also Section |5]. 
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For instance for the reconstruction problem of Example 1 with noise-free data and 
machine precision e = 5 = 1CT 16 , formula (pl|) implies r C [10 _6 ,10~ 8 ]. This coincides 
very well with numerical experiments. 

If the noise level 5 is not known, it has to be estimated. This difficult problem will 
not be discussed here. The reader is referred to p9f for more details. 



Although we have arrived now at an implementable algorithm for the nonuniform 
sampling problem, the disadvantages of the approach described in the previous section 
are obvious. In general the matrix R n does not have any particular structure, thus the 
computational costs for the singular value decomposition are 0(n 3 ) which is prohibitive 
large in many applications. It is definitely not a good approach to transform a well- 
posed infinite-dimensional problem into an ill-posed finite-dimensional problem for which 
a stable solution can only be computed by using a "heavy regularization machinery" . 

The methods in |42| , p| |40| , |33| , coincide with or are essentially equivalent to the 
truncated frame approach, therefore they suffer from the same instability problems and 
the same numerical inefficiency. 

2.2 CG and regularization of the truncated frame method 

As mentioned above one way to stabilize the solution of R n c^ = is a truncated sin- 
gular value decomposition, where the truncation level serves as regularization parameter. 
For large n the costs of the singular value decomposition become prohibitive for practical 
purposes. 

We propose the conjugate gradient method to solve R n S n ^ = b^ n \ It is in general 



much more efficient than a TSVD (or Tikhonov regularization as suggested in |^0|), and 
at the same time it can also be used cts db rc ffularization method. 

The standard error analysis for CG cannot be used in our case, since the matrix is 
ill-conditioned. Rather we have to resort to the error analysis developed in ||28| , |22fl . 

When solving a linear system Ax = y by CG for noisy data y s following happens. The 
iterates Xk of CG may diverge for k — > oo, however the error propagation remains limited 
in the beginning of the iteration. The quality of the approximation therefore depends on 
how many iterative steps can be performed until the iterates turn to diverge. The idea 
is now to stop the iteration at about the point where divergence sets in. In other words 
the iterations count is the regularization parameter which remains to be controlled by an 
appropriate stopping rule E7I, 



In our case assume 

U&M) _ < S\\b^\\, where &$ n,a) denotes a noisy sample. We 
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terminate the CG iterations when the iterates (c <n ' <5 ))/ c satisfy for the first time [T2 

0n) _ ( c W)) fc || < T S\\b^\\ (22) 

for some fixed r > 1. 

It should be noted that one can construct "academic" examples where this stopping 
rule does not prevent CG from diverging, see ||22|| , "most of the time" however it gives 



satisfactory results. We refer the reader to |27], |22| for a detailed discussion of various 
stopping criteria. 

There is a variety of reasons, besides the ones we have already mentioned, that make 
the conjugate gradient method and the nonuniform sampling problem a "perfect couple". 
See Sections |3|, [O], and £|7^ for more details. 

By combining the truncated frame approach with the conjugate gradient method (with 
appropriate stopping rule) we finally arrive at a reconstruction method that is of some 
practical relevance. However the only existing method at the moment that can handle 
large scale reconstruction problems seems to be the one proposed in the next section. 



3 Trigonometric polynomials and efficient signal re- 
construction 

In the previous section we have seen that the naive finite-dimensional approach via trun- 
cated frames is not satisfactory, it already leads to severe stability problems in the ideal 
case of regular oversampling. In this section we propose a different finite-dimensional 
model, which resembles much better the structural properties of the sampling problem, 
as can be seen below. 

The idea is simple. In practice only a finite number of samples {f{tj)} r j = i is given, 
where without loss of generality we assume — M < t\ < • • ■ < t r < M (otherwise we can 
always re-normalize the data). Since no data of / are available from outside this region we 
focus on a local approximation of / on [— M, M]. We extend the sampling set periodically 
across the boundaries, and identify this interval with the (properly normalized) torus T. 
To avoid technical problems at the boundaries in the sequel we will choose the interval 
somewhat larger and consider either [-M - 1/2, M + 1 /2] or [-N, N] with N = M + ^j. 
For theoretical considerations the choice [— M — 1/2, M + 1/2] is more convenient. 

Since the dual group of the torus T is Z, periodic band-limited functions on T reduce to 
trigonometric polynomials (of course technically / does then no longer belong to B since 
it is no longer in X 2 (R)). This suggests to use trigonometric polynomials as a realistic 
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finite-dimensional model for a numerical solution of the nonuniform sampling problem. 
We consider the space Pm of trigonometric polynomials of degree M of the form 

M 

p(t) = (2M + l)" 1 a k e 27Tikt/{2M+1) . (23) 

k=-M 

The norm of p G Pm is 



N M 

\\p\\ 2 = / \p(t)\ 2 dt= 



\a-k 



2 



-iV 



k=-M 



Since the distributional Fourier transform of p is p = (2M + 1) 1 J2k=-M a fc^/(2M+i) we 
have suppp C {k/(2M + l),\k\ < M} C [-1/2,1/2]. Hence P M is indeed a natural 
finite-dimensional model for B. 

In general the f(tj) are not the samples of a trigonometric polynomial in Pm, moreover 
the samples are usually perturbed by noise, hence we may not find a p G Pm such that 
p(tj) = bj = f(tj). We therefore consider the least squares problem 

r 

Here the lOj > are user-defined weights, which can be chosen for instance to compensate 
for irregularities in the sampling geometry WA . 



By increasing M so that r < 2M + 1 we can certainly find a trigonometric polynomial 
that interpolates the given data exactly. However in the presence of noise, such a solution 
is usually rough and highly oscillating and may poorly resemble the original signal. We 
will discuss the question of the optimal choice of M if the original bandwidth is not known 
and in presence of noisy data in Section |4.2| . 

The following theorem provides an efficient numerical reconstruction algorithm. It is 
also the key for the analysis of the relation between the finite-dimensional approximation 
in Pm and the solution of the original infinite-dimensional sampling problem in B. 



Theorem 3.1 (and Algorithm) / [J^ , [7^|/ Given the sampling points —M < t\ < . . . , t r < 

M , samples {bj} r j =1 , positive weights {vjj}j =1 with 2M + 1 < r. 

Step 1: Compute the (2M + 1) x (2M + 1) Toeplitz matrix Tm with entries 



r 

(T M )k,i = £ w .e-=M fc -%/(™+i) f 0r | fc | ? |/| < M (25) 

3=1 
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and y M E C< 2M+1 > by 



1 r 

(y M )k = ; y & 7 - Wj - e - 2,rifct ' /(2M+1) /or Ifcl < M . (26) 

V ; V2M + 1 ^ J J IU K ' 

Step 2: Solve the system 

T M aM = yn • (27) 



Step 3: Then the polynomial pu £ Pm that solves (|24| ) zs given by 

M 

PM(t) = -=^= £ M^* /(2¥+1) • (28) 
V2M + 1 fc= _ M 



Numerical Implementation of Theorem/ Algorithm |3.1j 



Step 1: The entries of T M and yu of equations ( 551 ) and (EH) can be computed in 
(9(MlogM + rlog(l/e)) operations (where e is the required accuracy) using Beylkin's 
unequally spaced FFT algorithm 0]. 

Step 2: We solve TmOm = yM by the conjugate gradient (CG) algorithm [^8[ . The 
matrix- vector multiplication in each iteration of CG can be carried out in O(MlogM) 
operations via FFT ||. Thus the solution of ( p7|) takes 0{kM log M) operations, where 
k is the number of iterations. 

Step 3: Usually the signal is reconstructed on regularly space nodes {-Uj}^. In this case 
PM{ui) in (|28|) can be computed by FFT. For non- uniformly spaced nodes U{ we can again 
resort to Beylkin's USFFT algorithm. 

There exists a large number of fast algorithms for the solution of Toeplitz systems. 
Probably the most efficient algorithm in our case is CG. We have already mentioned that 
the Toeplitz system ( |2"T| ) can be solved in 0(kM log M) via CG. The number of iterations 
k depends essentially on the clustering of the eigenvalues of Tm, cf. ||. It follows from 
equation (|3~I| ) below and perturbation theory [IU that, if the sampling points stem from 



a perturbed regular sampling set, the eigenvalues of Tm will be clustered around /3, where 
(3 is the oversampling rate. In such cases we can expect a very fast rate of convergence. 
The simple frame iteration |2B, |T| is not able to take advantage of such a situation. 



For the analysis of the relation between the solution pm of Theorem |3.1| and the solution 
/ of the original infinite-dimensional problem we follow Grochenig |2Hj. Assume that 
the samples {f{tj)}j^z of / G B are given. For the finite-dimensional approximation we 
consider only those samples f(tj) for which tj is contained in the interval [— M — |, M + 1] 
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and compute the least squares approximation pu with degree M and period 2M + 1 as 
in Theorem |3.1| . It is shown in |2(| that if <j(Tm) Q [a, /3] for all M with a > then 

Km / \f(t)- PM (t)\ 2 dt = 0, (29) 

[~M,M] 

and also YmvpM (t) = f(t) uniformly on compact sets. 

Under the Nyquist condition sup(i, + i — tj) := 7 < 1 and using weights Wj = (tj+i — 
tj-i)/2 Grochenig has shown that 

ct(T m ) C [(l- 7 ) 2 ,6], (30) 



independently of M, see ||20|| . These results validate the usage of trigonometric polyno- 
mials as finite-dimensional model for nonuniform sampling. 

Example 1 — reconsidered: Recall that in Example 1 of Section |2| we have considered 
the reconstruction of a regularly oversampled signal / G B. What does the reconstruction 
method of Theorem |3. 1| yield in this case? Let us check the entries of the matrix Tm when 
we take only those samples in the interval [— n, n}. The period of the polynomial becomes 
2N with N = n + -Aj- where r is the number of given samples. Then 

r nm 

(TmU = ^ E e 2 ^^ = J2 e 2 ^-<>™ = m6 k>l (31) 

j=l j=—nm 

for k, I = —M, . . . , M, where 8^1 is Kronecker's symbol with the usual meaning 8j.j = 1 if 
k = I and else. Hence we get 

T M = ml , 

where I is the identity matrix on C 2M+1 , thus Tm resembles the structure of the infinite- 
dimensional frame operator S in this case (including exact approximation of the frame 
bounds). Recall that the truncated frame approach leads to an "artificial" ill-posed prob- 
lem even in such a simple situation. 

The advantages of the trigonometric polynomial approach compared to the truncated 
frame approach are manifold. In the one case we have to deal with an ill-posed problem 
which has no specific structure, hence its solution is numerically very expensive. In the 
other case we have to solve a problem with rich mathematical structure, whose stability 
depends only on the sampling density, a situation that resembles the original infinite- 
dimensional sampling problem. 
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In principle the coefficients a M = {{dM)k\h=-M °f the polynomial pm that mini- 
mizes (|24]) could also be computed by directly solving the Vandermonde type system 

WVa M = Wb , (32) 

where V jjk = _l^ e -»i/(^+ 1 ) f or j = 1, . . . , r, k — -M, ...,M and W is a diagonal 



matrix with entries Wjj = y/w], cf. ||31|| . Several algorithms are known for a relatively 
efficient solution of Vandermonde systems [||, [3T|. However this is one of the rare cases, 
where, instead of directly solving (|32|), it is advisable to explicitly establish the system of 
normal equations 

Tm&m = Vm , (33) 

where T = V*W 2 V and y = V*W 2 b. 

The advantages of considering the system Tmcim = Vm instead of the Vandermonde 
system ([32] ) are manifold: 

• The matrix Tm plays a key role in the analysis of the relation of the solution of (|24|) 
and the solution of the infinite-dimensional sampling problem (||) , see ( ^9|) and (130|) 
above. 

• Tm is of size (2M + 1) x (2M+ 1), independently of the number of sampling points. 
Moreover, since (T M ) k j = Y^j=i Wje 2m( - k ~ l,)t: > ', it is of Toeplitz type. These facts give 
rise to fast and robust reconstruction algorithms. 

• The resulting reconstruction algorithms can be easily generalized to higher dimen- 



sions, see Section 3.1. Such a generalization to higher dimensions seems not to 



be straightforward for fast solvers of Vandermonde systems such as the algorithm 



proposed in |3T| 



We point out that other finite-dimensional approaches are proposed in |L6], [7|. These 
approaches may provide interesting alternatives in the few cases where the algorithm 
outlined in Section |] does not lead to good results. These cases occur when only a few 
samples of the signal / are given in an interval [a, b] say, and at the same time we have 
\f(a) - f(b)\ > and \f'(a) - f'(b)\ > 0, i.e., if / is "strongly non-periodic" on [a, b). 
However the computational complexity of the methods in |16], is significantly larger. 

3.1 Multi-dimensional nonuniform sampling 

The approach presented above can be easily generalized to higher dimensions by a dili- 
gent book-keeping of the notation. We consider the space of c/-dimensional trigonometric 
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polynomials P M as finite-dimensional model for B . For given samples f(tj) of / & B , 
where tj G M. d , we compute the least squares approximation pm similar to Theorem [TT 
by solving the corresponding system of equations Tmclm = Vm- 

In 2-D for instance the matrix Tm becomes a block Toeplitz matrix with Toeplitz 
blocks |37j . For a fast computation of the entries of T we can again make use of Beylkin's 



USFFT algorithm And similar to 1-D, multiplication of a vector by Tm can be carried 
out by 2-D FFT. 

Also the relation between the finite-dimensional approximation in P M and the infinite- 
dimensional solution in B d is similar as in 1-D. The only mathematical difficulty is to 
give conditions under which the matrix Tm is invertible. Since the fundamental theorem 
of algebra does not hold in dimensions larger than one, the condition (2M + l) d < r is 
necessary but no longer sufficient for the invertibility of Tm- Sufficient conditions for the 



invertibility, depending on the sampling density, are presented in [21 



4 Bandwidth estimation and regularization 

In this section we discuss several numerical aspects of nonuniform sampling that are very 
important from a practical viewpoint, however only few answers to these problems can 
be found in the literature. 



4.1 A multilevel signal reconstruction algorithm 

In almost all theoretical results and numerical algorithms for reconstructing a band-limited 
signal from nonuniform samples it is assumed that the bandwidth is known a priori. This 
information however is often not available in practice. 

A good choice of the bandwidth for the reconstruction algorithm becomes crucial in 
case of noisy data. It is intuitively clear that choosing a too large bandwidth leads to 
over-fit of the noise in the data, while a too small bandwidth yields a smooth solution but 
also to under-fit of the data. And of course we want to avoid the determination of the 
"correct" Q by trial-and-error methods. Hence the problem is to design a method that 
can reconstruct a signal from non-uniformly spaced, noisy samples without requiring a 
priori information about the bandwidth of the signal. 

The multilevel approach derived in provides an answer to this problem. The ap- 



proach applies to an infinite-dimensional as well as to a finite-dimensional setting. We 
describe the method directly for the trigonometric polynomial model, where the determi- 
nation of the bandwidth Q translates into the determination of the polynomial degree M 
of the reconstruction. The idea of the multilevel algorithm is as follows. 
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Let the noisy samples {£>*}£ =1 = {f s (tj)} r j =1 of / G B be given with Y^j=i \f{tj) " 
{tj)\ 2 < 5 2 \\b s \\ 2 and let Qm denote the orthogonal projection from B into Pm- We 



start with initial degree M = 1 and run Algorithm [3.1| until the iterates po,fc satisfy for 
the first time the inner stopping criterion 

r 

£ \ Pl , k (tj) - bf < 2r(5\\b 5 \\ + \\Q f - f\\)\\b s \\ 
j=i 

for some fixed r > 1. Denote this approximation (at iteration k*) by Pi,fc„. If Pi.fc, satisfies 
the cmter stopping criterion 

r 

^2\pi, k (tj)-b S j \ 2 <2r8\\b s \\ 2 (34) 
we take p\ t k t as final approximation. Otherwise we proceed to the next level M = 2 and 



run Algorithm [D] again, using as initial approximation by setting p 2 ,o = 



At level M = N the inner level-dependent stopping criterion becomes 

r 

\PN,k(tj) ~ b 5 j\ 2 < 2r(5\\b 5 \\ + \\Q N f - f\\)\\b s \\, (35) 

3=1 

while the outer stopping criterion does not change since it is level- independent. 

Stopping rule (|35|) guarantees that the iterates of CG do not diverge. It also ensures 
that CG does not iterate too long at a certain level, since if M is too small further 
iterations at this level will not lead to a significant improvement. Therefore we switch 
to the next level. The outer stopping criterion (34) controls over-fit and under-fit of the 
data, since in presence of noisy data is does not make sense to ask for a solution pu that 
satisfies Y7j=i Ipm(^) - = 0. 

Since the original signal / is not known, the expression ||/ — Qwf\\ in fl3"5]) cannot be 



computed. In |34| the reader can find an approach to estimate ||/ — QnJW recursively. 



4.2 Solution of ill-conditioned sampling problems 

A variety of conditions on the sampling points {tj}j £ ^ are known under which the set 
{sinc(- —tj)}j e z is a frame for B, which in turn implies (at least theoretically) perfect 
reconstruction of a signal / from its samples f(tj). This does however not guarantee 
a stable reconstruction from a numerical viewpoint, since the ratio of the frame bounds 
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B/A can still be extremely large and therefore the frame operator 5* can be ill-conditioned. 
This may happen for instance if 7 in ( p0[) goes to 1, in which case cond(T) may become 
large. The sampling problem may also become numerically unstable or even ill-posed, 
if the sampling set has large gaps, which is very common in astronomy and geophysics. 
Note that in this case the instability of the system Tm&m — Vm does not result from an 
inadequate discretization of the infinite-dimensional problem. 

There exists a large number of (circulant) Toeplitz preconditioners that could be ap- 
plied to the system Tmclm = Vm, however it turns out that they do not improve the 
stability of the problem in this case. The reason lies in the distribution of the eigenvalues 
of Tm, as we will see below. 

Following |38[1 , we call two sequences of real numbers {\^}^ =1 and {^ n '}k=i equally 
distributed, if 



lim 

n—>oo n 



n 

-E[ F ( A i n) )- F ^ n) )] = ° ( 36 ) 



k=l 



for any continuous function F with compact support^. 

Let C be a (n x n) circulant matrix with first column (co, . . . , c n _i), we write C = 
circ(co, . . . , c n _i). The eigenvalues of C are distributed as = Ym=o qe 27nfc///n . Ob- 
serve that the Toeplitz matrix A n with first column (a , at, ... , a n ) can be embedded in 
the circulant matrix 

C n = circ(a , a 1: . . . , a n , a n , . . . , ai) . (37) 



Thms 4.1 and 4.2 in ]38| state that the eigenvalues of A n and C n are equally distributed 
as f(x) where 



00 



fix) = a ^ 2mkx ■ (38) 

k=— 00 

The partial sum of the series ( |3"8"D is 

n 

f n (x)= a k e 2mkx . (39) 



k=—n 



To understand the clustering behavior of the eigenvalues of Tm in case of sampling 
sets with large gaps, we consider a sampling set in [— M, M), that consists of one large 



x In H.Weyl's definition a[™' and v£"' are required to belong to a common interval. 
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block of samples and one large gap, i.e., tj = j— for j = —mM, . . . mM for m, L E N. 
(Recall that we identify the interval with the torus). Then the entries Zk of the Toeplitz 
matrix T M of fl25|) (with uij = 1) are 

mM 



I Y e -2^i/(2M+l) k = Q 2M. 

+ 1 ^ 



Zk 2M , 

j=—mM 

To investigate the clustering behavior of the eigenvalues of Tm for M— > oo, we embed Tj 



in a circulant matrix Cm as in (|37|) . Then (p9T) becomes 

Ul-)- E E e-«(«'+l>-«(-+l>->l (40) 

^ ' l=-mM j=-mM 

whence f m M — > l[-i/(2L),i/(2L)] for oo, where l[_ a](1 ](x) = 1, if —a < x < a and else. 

Thus the eigenvalues of T M are asymptotically clustered around zero and one. For 
general nonuniform sampling sets with large gaps the clustering at 1 will disappear, but of 
course the spectral cluster at will remain. In this case it is known that the preconditioned 



problem will still have a spectral cluster at the origin [|43| and preconditioning will not be 
efficient. 

Fortunately there are other possibilities to obtain a stabilized solution of Tmo>m = Um- 
The condition number of Tm essentially depends on the ratio of the maximal gap in 
the sampling set to the Nyquist rate, which in turn depends on the bandwidth of the 
signal. We can improve the stability of the system by adapting the degree M of the 
approximation accordingly. Thus the parameter M serves as a regularization parameter 
that balances stability and accuracy of the solution. This technique can be seen as a 
specific realization of regularization by projection, see Chapter 3 in f[2| . In addition, as 
described in Section [4.2| , we can utilize CG as regularization method for the solution of 
the Toeplitz system in order to balance approximation error and propagated error. The 



multilevel method introduced in Section 4.1 combines both features. By optimizing the 



level (bandwidth) and the number of iterations in each level it provides an efficient and 
robust regularization technique for ill-conditioned sampling problems. See Section |^ for 
numerical examples. 



5 Applications 

We present two numerical examples to demonstrate the performance of the described 
methods. The first one concerns a 1-D reconstruction problem arising in spectroscopy. In 
the second example we approximate the Earth's magnetic field from noisy scattered data. 
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5.1 An example from spectroscopy 

The original spectroscopy signal / is known at 1024 regularly spaced points tj. This 
discrete sampling sequence will play the role of the original continuous signal. To simulate 
the situation of a typical experiment in spectroscopy we consider only 107 randomly chosen 
sampling values of the given sampling set. Furthermore we add noise to the samples with 
noise level (normalized by division by Ylk=i l/(^i)| 2 ) of 5 = 0.1. Since the samples are 
contaminated by noise, we cannot expect to recover the (discrete) signal / completely. 
The bandwidth is approximately Q = 5 which translates into a polynomial degree of 
M ~ 30. Note that in general fl and (hence M) may not be available. We will also 
consider this situation, but in the first experiments we assume that we know Q. The 
error between the original signal / and an approximation /„ is measured by computing 

ll/-/n|| 2 /ll/H 2 - 

First we apply the truncated frame method with regularized SVD as described in 



Section |2|. We choose the truncation level for the SVD via formula (EJ]). This is the 
optimal truncation level in this case, providing an approximation with least squares error 
0.0944. Figure |l](a) shows the reconstructed signal together with the original signal and 
the noisy samples. Without regularization we get a much worse "reconstruction" (which 
is not displayed). 



We apply CG to the truncated frame method, as proposed in Section ^2] with stop- 
ping criterion (53) (for r = 1). The algorithm terminates already after 3 iterations. 
The reconstruction error is with 0.1097 slightly higher than for truncated SVD (see also 
Figure [j](b)), but the computational effort is much smaller. 

Also Algorithm [34] (with M = 30) terminates after 3 iterations. The reconstruction 
is shown in Figure 0(c), the least squares error (0.0876) is slightly smaller than for the 
truncated frame method, the computational effort is significantly smaller. 

We also simulate the situation where the bandwidth is not known a priori and demon- 
strate the importance of a good estimate of the bandwidth. We apply Algorithm |3.1| using 
a too small degree (M =11) and a too high degree (M = 40). (We get qualitatively the 
same results using the truncated frame method when using a too small or too large band- 
width). The approximations are shown in Figs. |l|(d) and (e), The approximation errors 



are 0.4648 and 0.2805, respectively. Now we apply the multilevel algorithm of Section [O 
which does not require any initial choice of the degree M. The algorithm terminates at 
"level" M = 22, the approximation is displayed in Fig. |l|(f), the error is 0.0959, thus 
within the error bound 5, as desired. Hence without requiring explicit information about 
the bandwidth, we are able to obtain the same accuracy as for the methods above. 
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(a) Truncated frame method (b) Truncated frame method 
with TSVD, error=0.0944. with CG, error=0.1097. 





* nois sam les 
, original signal 


V 






(c) Algorithm 3.1 with "correct" 
bandwidth, error=0.0876 



(d) Using a too small band- 
width, error=0.4645. 









1 \ 


1 ; 




(e) Using a too large bandwidth, 
error = 0.2412. 



(f) Multilevel algorithm, er- 
ror=0.0959. 



Figure 1: Example from spectroscopy - comparison of reconstruction methods. 
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5.2 Approximation of geophysical potential fields 

Exploration geophysics relies on surveys of the Earth's magnetic field for the detection 
of anomalies which reveal underlying geological features. Geophysical potential field-data 
are generally observed at scattered sampling points. Geoscientists, used to looking at 
their measurements on maps or profiles and aiming at further processing, therefore need 
a representation of the originally irregularly spaced data at a regular grid. 

The reconstruction of a 2-D signal from its scattered data is thus one of the first and 
crucial steps in geophysical data analysis, and a number of practical constraints such 
as measurement errors and the huge amount of data make the development of reliable 
reconstruction methods a difficult task. 

It is known that the Fourier transform of a geophysical potential field / has decay 
| f(u) | = (9(e~' w '). This rapid decay implies that / can be very well approximated by 
band-limited functions |3(J . Since in general we may not know the (essential) bandwidth 



of /, we can use the multilevel algorithm proposed in Section [4J] to reconstruct /. 

The multilevel algorithm also takes care of following problem. Geophysical sampling 
sets are often highly anisotropic and large gaps in the sampling geometry are very common. 
The large gaps in the sampling set can make the reconstruction problem ill-conditioned or 
even ill-posed. As outlined in Section ^72| the multilevel algorithm iteratively determines 
the optimal bandwidth that balances the stability and accuracy of the solution. 



Figure ^(a) shows a synthetic gravitational anomaly /. The spectrum of / decays 
exponentially, thus the anomaly can be well represented by a band-limited function, using 
a "cut-off-level" of |/(u;)| < 0.01 for the essential bandwidth of /. 

We have sampled the signal at 1000 points (v,j, Vj) and added 5% random noise to the 
sampling values f(v,j,Vj). The sampling geometry - shown in Figure [5.2| as black dots 
- exhibits several features one encounters frequently in exploration geophysics ||30|| . The 
essential bandwidth of / would imply to choose a polynomial degree of M = 12 (i.e., 
(2M + l) 2 = 625 spectral coefficients). With this choice of M the corresponding block 
Toeplitz matrix Tm would become ill-conditioned, making the reconstruction problem 
unstable. As mentioned above, in practice we usually do not know the essential bandwidth 
of /. Hence we will not make use of this knowledge in order to approximate /. 

We apply the multilevel method to reconstruct the signal, using only the sampling 
points {(uj,Vj)}, the samples {f 5 (uj, v j)} and the noise level 5 = 0.05 as a priori infor- 
mation. The algorithm terminates at level M — 7. The reconstruction is displayed in 
Figure |5.2| (c), the error between the true signal and the approximation is shown in Fig- 
ure |5.2|(d). The reconstruction error is 0.0517 (or 0.193 mGal), thus of the same order as 



the data error, as desired. 
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easting [m] easting [m] 



(a) Contour map of synthetic grav- (b) Sampling set and synthetic 

ity anomaly, gravity is in mGal. gravity anomaly. 




(c) Approximation by multi-level al- (d) Error between approximation 

gorithm and actual anomaly. 

Figure 2: Approximation of synthetic gravity anomaly from 1000 non-uniformly spaced 
noisy samples by the multilevel algorithm of Section [4.1| . The algorithm iteratively deter- 
mines the optimal bandwidth (i.e. level) for the approximation. 
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